#####################################################################################
# Alliance and Public Preference for Nuclear Forbearance: Evidence from South Korea 
#####################################################################################


rm(list=ls())

library("foreign")


###########################################################################
########### Table 1. The Dampening Effect of Three Non-proliferation Tools 
###########################################################################

##load the public data 


datap<-read.csv("public.csv")
data3p<-subset(datap, datap$age==1|datap$age==2|datap$age==3) #20-30s block
data4p<-subset(datap, datap$age==4|datap$age==5|datap$age==6) #40-50s block


###Bootstrp function 

B<-10000

bootTreat<- function(dat, label="Treatment"){
  bootResults <- matrix(NA, nrow=B, ncol=3)
  for (i in seq_len(B)){
    resample <- sample(1:nrow(dat),nrow(dat),replace=T)
    temp <- dat[resample,]
    A <- mean(temp$support[which(temp$declaration==1)], na.rm=TRUE)
    B <- mean(temp$support[which(temp$deployment==1)], na.rm=TRUE)
    C <- mean(temp$support[which(temp$threat==1)], na.rm=TRUE)
    D <- mean(temp$support[which(temp$control==1)], na.rm=TRUE)
    bootResults[i,1] <- (A-D)
    bootResults[i,2] <- (B-D)
    bootResults[i,3] <- (C-D)
    drop(list())
  }
  return(list(model=label, boot=bootResults, Declare=mean(bootResults[,1], na.rm=TRUE), Deploy=mean(bootResults[,2], na.rm=TRUE), Threat=mean(bootResults[,3], na.rm=TRUE)))
}


### full public sample 

set.seed(1000)
full.p <- bootTreat(datap)


mean(datap$support[datap$control==1])
mean(datap$support[datap$declaration==1])-mean(datap$support[datap$control==1])
mean(datap$support[datap$deploy==1])-mean(datap$support[datap$control==1])
mean(datap$support[datap$threat==1])-mean(datap$support[datap$control==1])

length(which(full.p$boot[,1] >0))/B #p=0.3261: declaration
length(which(full.p$boot[,2] >0))/B #p=0.1804: deployment
length(which(full.p$boot[,3] >0 ))/B #p=0.4077: threat

##simple t-test give you the same results

t.test(datap$support[datap$declaration==1], datap$support[datap$control==1], alternative=c("less"))
#p=0.3222
t.test(datap$support[datap$deployment==1], datap$support[datap$control==1], alternative=c("less"))
#p=0.1773
t.test(datap$support[datap$threat==1], data3p$support[datap$control==1], alternative=c("less"))
#p=0.4541


###20-30s sample
set.seed(1000)
full.3p <- bootTreat(data3p)


mean(data3p$support[data3p$control==1])
mean(data3p$support[data3p$declaration==1])-mean(data3p$support[data3p$control==1])
mean(data3p$support[data3p$deploy==1])-mean(data3p$support[data3p$control==1])
mean(data3p$support[data3p$threat==1])-mean(data3p$support[data3p$control==1])


length(which(full.3p$boot[,1]>0))/B #p=0.0517
length(which(full.3p$boot[,2]>0))/B #p=0.034
length(which(full.3p$boot[,3]>0 ))/B #p=0.0138




t.test(data3p$support[data3p$declaration==1], data3p$support[data3p$control==1], alternative=c("less"))
#p=0.05152
t.test(data3p$support[data3p$deployment==1], data3p$support[data3p$control==1], alternative=c("less"))
#p=0.03619
t.test(data3p$support[data3p$threat==1], data3p$support[data3p$control==1], alternative=c("less"))
#p=0.01644


####40-50s sample

set.seed(1000)
full.4p <- bootTreat(data4p)


mean(data4p$support[data4p$control==1])
mean(data4p$support[data4p$declaration==1])-mean(data4p$support[data4p$control==1])
mean(data4p$support[data4p$deploy==1])-mean(data4p$support[data4p$control==1])
mean(data4p$support[data4p$threat==1])-mean(data4p$support[data4p$control==1])



length(which(full.4p$boot[,1]<0))/B #p=0.1568
length(which(full.4p$boot[,2]<0))/B #p=0.311
length(which(full.4p$boot[,3]<0 ))/B #p=0.0343


## T-test 
t.test(data4p$support[data4p$declaration==1], data4p$support[data4p$control==1], alternative=c("greater"))
#p=0.1642
t.test(data4p$support[data4p$deployment==1], data4p$support[data4p$control==1], alternative=c("greater"))
#p=0.3119
t.test(data4p$support[data4p$threat==1], data4p$support[data4p$control==1], alternative=c("greater"))
#p=0.0366



#######################################################################################################
######### Table 2. The Dampening Effect of Three Non-proliferation Tools in the Survey on Government Officials
########################################################################################################

datae<-read.csv("elite.csv")


mean(datae$support[datae$control==1])
mean(datae$support[datae$declaration==1])-mean(datae$support[datae$control==1])
mean(datae$support[datae$deploy==1])-mean(datae$support[datae$control==1])
mean(datae$support[datae$threat==1])-mean(datae$support[datae$control==1])



set.seed(1000)
full.e <-bootTreat(datae)

length(which(full.e$boot[,1] >0))/B #p=0.0178
length(which(full.e$boot[,2] >0))/B #p=0.0146
length(which(full.e$boot[,3] >0 ))/B #p=0.1222




## T-test 
t.test(datae$support[datae$declaration==1], datae$support[datae$control==1], alternative=c("less"))
#p=0.02204
t.test(datae$support[datae$deployment==1], datae$support[datae$control==1], alternative=c("less"))
#p=0.0211
t.test(datae$support[datae$threat==1], datae$support[datae$control==1], alternative=c("less"))
#p=0.1424




###################################################################################################
#### To Create Figure 1 and 2 calculate the difference in means and confidence intervals (Public Survey)
####################################################################################################


### Bootstrap Functions

B<-10000

bootTreat2 <- function(dat, label="Treatment"){
  bootResults <- matrix(NA, nrow=B, ncol=2)
  for (i in seq_len(B)){
    resample <- sample(1:nrow(dat),nrow(dat),replace=T)
    temp <- dat[resample,]
    A <- mean(temp$support[which(temp$declaration==1)], na.rm=TRUE)
    B <- mean(temp$support[which(temp$deployment==1)], na.rm=TRUE)
    C <- mean(temp$support[which(temp$threat==1)], na.rm=TRUE)
    D <- mean(temp$support[which(temp$control==1)], na.rm=TRUE)
    bootResults[i,1] <- (B-A)
    bootResults[i,2] <- (C-A)
    drop(list())
  }
  return(list(model=label, boot=bootResults, Deploy=mean(bootResults[,1], na.rm=TRUE), Threat=mean(bootResults[,2], na.rm=TRUE)))
}


bootCredit2 <- function(dat, label="Treatment"){
  bootResults <- matrix(NA, nrow=B, ncol=2)
  for (i in seq_len(B)){
    resample <- sample(1:nrow(dat),nrow(dat),replace=T)
    temp <- dat[resample,]
    A <- mean(temp$credibility[which(temp$declaration==1)], na.rm=TRUE)
    B <- mean(temp$credibility[which(temp$deployment==1)], na.rm=TRUE)
    C <- mean(temp$credibility[which(temp$threat==1)], na.rm=TRUE)
    bootResults[i,1] <- (B-A)
    bootResults[i,2] <- (C-A)
    drop(list())
  }
  return(list(model=label, boot=bootResults, Deploy=mean(bootResults[,1], na.rm=TRUE), Threat=mean(bootResults[,2], na.rm=TRUE)))
}



####20s-30s block 

##1) support for nuclearization 


mean(data3p$support[data3p$deploy==1])-mean(data3p$support[data3p$declaration==1]) #0.013
mean(data3p$support[data3p$threat==1])-mean(data3p$support[data3p$declaration==1]) #0.04

##calculate confidence intervals 
set.seed(1000)
full.3pr <- bootTreat2(data3p)

p3d<-round(c(full.3pr$Deploy, quantile(full.3pr$boot[,1], c(0.025, 0.975))), digits=2)
p3d # (-0.17, 0.14)

p3t<-round(c(full.3pr$Threat, quantile(full.3pr$boot[,2], c(0.025, 0.975))), digits=2)
p3t # (-0.20, 0.11)

p3d90<-round(c(full.3pr$Deploy, quantile(full.3pr$boot[,1], c(0.05, 0.95))), digits=2)
p3d90 # (-0.14, 0.12)

p3t90<-round(c(full.3pr$Threat, quantile(full.3pr$boot[,2], c(0.05, 0.95))), digits=2)
p3t90 #(-0.17, 0.09)


##2) credibility of extended deterrence 
mean(data3p$credibility[data3p$declaration==1]) #0.213

mean(data3p$credibility[data3p$deploy==1])-mean(data3p$credibility[data3p$declaration==1]) #0.106
mean(data3p$credibility[data3p$threat==1])-mean(data3p$credibility[data3p$declaration==1]) #0.133


 #calculate confidence intervals 
set.seed(1000)
full.3pc<-bootCredit2(data3p)


c3d<-round(c(full.3pc$Deploy, quantile(full.3pc$boot[,1], c(0.025, 0.975))), digits=2)
c3d #0.11 (-0.04, 0.24)

c3t<-round(c(full.3pc$Threat, quantile(full.3pc$boot[,2], c(0.025, 0.975))), digits=2)
c3t #0.13 (-0.01, 0.28)

c3d90<-round(c(full.3pc$Deploy, quantile(full.3pc$boot[,1], c(0.05, 0.95))), digits=2)
c3d90 #0.11 (-0.01, 0.22)

c3t90<-round(c(full.3pc$Threat, quantile(full.3pc$boot[,2], c(0.05, 0.95))), digits=2)
c3t90 #0.13(0.01, 0.25)



####40s-50s Block 



##1) support for nuclearization 


mean(data4p$support[data4p$deploy==1])-mean(data4p$support[data4p$declaration==1]) #-0.04
mean(data4p$support[data4p$threat==1])-mean(data4p$support[data4p$declaration==1]) #0.066


###calculate confidence intervals
set.seed(1000)
full.4pr <- bootTreat2(data4p)

p4d<-round(c(full.4pr$Deploy, quantile(full.4pr$boot[,1], c(0.025, 0.975))), digits=2)
p4d# (-0.20, 0.12)

p4t<-round(c(full.4pr$Threat, quantile(full.4pr$boot[,2], c(0.025, 0.975))), digits=2)
p4t# (-0.09, 0.23)

p4d90<-round(c(full.4pr$Deploy, quantile(full.4pr$boot[,1], c(0.05, 0.95))), digits=2)
p4d90 #(-0.17, 0.09)

p4t90<-round(c(full.4pr$Threat, quantile(full.4pr$boot[,2], c(0.05, 0.95))), digits=2)
p4t90 #(-0.07, 0.20)



##2) credibility of extended deterrence 
mean(data4p$credibility[data4p$declaration==1]) #0.28

mean(data4p$credibility[data4p$deploy==1])-mean(data4p$credibility[data4p$declaration==1]) #0.08
mean(data4p$credibility[data4p$threat==1])-mean(data4p$credibility[data4p$declaration==1]) #0.013


###calculate confidence intervals
set.seed(1000)
full.4pc<-bootCredit2(data4p)

c4d<-round(c(full.4pc$Deploy, quantile(full.4pc$boot[,1], c(0.025, 0.975))), digits=2)
c4d#(-0.07, 0.23)

c4t<-round(c(full.4pc$Threat, quantile(full.4pc$boot[,2], c(0.025, 0.975))), digits=2)
c4t #(-0.13, 0.16)

c4d90<-round(c(full.4pc$Deploy, quantile(full.4pc$boot[,1], c(0.05, 0.95))), digits=2)
c4d90 #(-0.05, 0.20)

c4t90<-round(c(full.4pc$Threat, quantile(full.4pc$boot[,2], c(0.05, 0.95))), digits=2)
c4t90 #(-0.11, 0.14)




### full sample 



##1) support for nuclearization 


mean(datap$support[datap$deploy==1])-mean(datap$support[datap$declaration==1]) ##-0.026
mean(datap$support[datap$threat==1])-mean(datap$support[datap$declaration==1]) ##0.013


##calculate confidence intervals
set.seed(1000)
full.pr <- bootTreat2(datap)


pd<-round(c(full.pr$Deploy, quantile(full.pr$boot[,1], c(0.025, 0.975))), digits=2)
pd# (-0.14,0.08)

pt<-round(c(full.pr$Threat, quantile(full.pr$boot[,2], c(0.025, 0.975))), digits=2)
pt# (-0.10, 0.13)

pd90<-round(c(full.pr$Deploy, quantile(full.pr$boot[,1], c(0.05, 0.95))), digits=2)
pd90# (-0.12, 0.07)

pt90<-round(c(full.pr$Threat, quantile(full.pr$boot[,2], c(0.05, 0.95))), digits=2)
pt90# (-0.08, 0.11)



##2) credibility of extended deterrence 
mean(datap$credibility[datap$declaration==1]) #0.246

mean(datap$credibility[datap$deploy==1])-mean(datap$credibility[datap$declaration==1]) ##0.093
mean(datap$credibility[datap$threat==1])-mean(datap$credibility[datap$declaration==1]) ##0.073


###calculate confidence intervals
set.seed(1000)
full.pc<-bootCredit2(datap)

cd<-round(c(full.pc$Deploy, quantile(full.pc$boot[,1], c(0.025, 0.975))), digits=2)
cd # (-0.01, 0.20)
ct<-round(c(full.pc$Threat, quantile(full.pc$boot[,2], c(0.025, 0.975))), digits=2)
ct #(-0.03, 0.18)
cd90<-round(c(full.pc$Deploy, quantile(full.pc$boot[,1], c(0.05, 0.95))), digits=2)
cd90 #(0.01, 0.18)
ct90<-round(c(full.pc$Threat, quantile(full.pc$boot[,2], c(0.05, 0.95))), digits=2)
ct90 # (-0.01, 0.16)




###combine results from the public survey into one dataframe



dataret<-as.data.frame(cbind(c("40/50 Block", "20/30 Block","Full Sample")))
dataret$dif<-cbind(c(p4t[1], p3t[1], pt[1] ))
dataret$low<-cbind(c(p4t[2], p3t[2], pt[2] ))
dataret$up<-cbind(c(p4t[3], p3t[3], pt[3] ))
dataret$low90<-cbind(c(p4t90[2], p3t90[2], pt90[2]))
dataret$up90<-cbind(c(p4t90[3], p3t90[3], pt90[3]))


dataret$dif<-dataret$dif*100
dataret$low<-dataret$low*100
dataret$up<-dataret$up*100
dataret$low90<-dataret$low90*100
dataret$up90<-dataret$up90*100



datared<-as.data.frame(cbind(c("40/50 Block", "20/30 Block","Full Sample" )))
datared$dif<-cbind(c(p4d[1], p3d[1], pd[1] ))
datared$low<-cbind(c(p4d[2], p3d[2], pd[2]))
datared$up<-cbind(c(p4d[3], p3d[3], pd[3]))
datared$low90<-cbind(c(p4d90[2], p3d90[2], pd90[2]))
datared$up90<-cbind(c(p4d90[3], p3d90[3], pd90[3]))


datared$dif<-datared$dif*100
datared$low<-datared$low*100
datared$up<-datared$up*100
datared$low90<-datared$low90*100
datared$up90<-datared$up90*100


datacet<-as.data.frame(cbind(c("40/50 Block", "20/30 Block","Full Sample")))
datacet$dif<-cbind(c(c4t[1], c3t[1], ct[1] ))
datacet$low<-cbind(c(c4t[2], c3t[2], ct[2]))
datacet$up<-cbind(c(c4t[3], c3t[3], ct[3]))
datacet$low90<-cbind(c(c4t90[2], c3t90[2], ct90[2]))
datacet$up90<-cbind(c(c4t90[3], c3t90[3], ct90[3]))

datacet$dif<-datacet$dif*100
datacet$low<-datacet$low*100
datacet$up<-datacet$up*100
datacet$low90<-datacet$low90*100
datacet$up90<-datacet$up90*100


dataced<-as.data.frame(cbind(c("40/50 Block", "20/30 Block","Full Sample")))
dataced$dif<-cbind(c(c4d[1], c3d[1], cd[1] ))
dataced$low<-cbind(c(c4d[2], c3d[2], cd[2]))
dataced$up<-cbind(c(c4d[3], c3d[3], cd[3]))
dataced$low90<-cbind(c(c4d90[2], c3d90[2], cd90[2]))
dataced$up90<-cbind(c(c4d90[3], c3d90[3], cd90[3]))

dataced$dif<-dataced$dif*100
dataced$low<-dataced$low*100
dataced$up<-dataced$up*100
dataced$low90<-dataced$low90*100
dataced$up90<-dataced$up90*100






###################################################################################
# Figure 1 The Relative EFfectivenss: Nuclear Foward-Deployment vs Declaratory Policy
####################################################################################



png("Figure1_FPA.png", width=10, height=6,  units = 'in', res = 300)

m<-rbind(c(1,2))
layout(m)
layout.show(2)

## Panel (a) Support for Nuclearization 
plot(NULL,                         
     xlim = c(-80, 80),                     
     ylim = c(.3, 4 + .3), 
     axes = F, xlab = NA, ylab = NA)      

abline(v = 0, lty = 2, col="grey") 

  
for (i in 1:length(datared$dif)) {                                          
  points(datared$dif[i], i, pch = 18, cex = 1.5)                              
  lines(c(datared$low[i], datared$up[i]), c(i, i))
  lines(c(datared$low90[i], datared$up90[i]),c(i,i), lwd=3)
  text(datared$dif[i], i, datared$V1[i], xpd = T, cex = 1, pos = 3)                    
}       

axis(side = 1)                                                                                    
mtext(side = 3, "(a) Support for Nuclearization", line = 0.9, cex=1.1)
box(col="grey")        


## Panel(b) Crediblity of Extended Deterrence

plot(NULL,                              
     xlim = c(-80, 80),                     
     ylim = c(.3, 4 + .3), 
     axes = F, xlab = NA, ylab = NA)       

abline(v = 0, lty = 2, col="grey") 
               
for (i in 1:length(dataced$dif)) {                                           
  points(dataced$dif[i], i, pch = 19, cex = 1.2)                              
  lines(c(dataced$low[i], dataced$up[i]), c(i, i))
  lines(c(dataced$low90[i], dataced$up90[i]),c(i,i), lwd=3)
  text(dataced$dif[i], i, dataced$V1[i], xpd = T, cex = 1, pos = 3) 
}

axis(side = 1)                                                                                         
mtext(side = 3, "(b) Credibility of Extended Deterrence ", line = 0.9, cex=1.1)

mtext(side = 1, "Changes in Support for Nuclearization \n(Deploy-Declare)(%)", line = 3, at=-220, cex=1) 
mtext(side = 1, "Changes in Credibility of Extended Deterrence \n(Deploy-Declare)(%)", line = 3, at=10, cex=1) 
box(col="grey")               


dev.off()




################################################################################
# Figure 2 The Relative EFfectivenss: a Conditional Threat of Punishment vs Declaratory Policy
#################################################################################




png("Figure2_FPA.png", width=10, height=6,  units = 'in', res = 300)


m<-rbind(c(1,2))
layout(m)

#Panel(a) Support for Nuclearization 


plot(NULL,                              
     xlim = c(-80, 80),                     
     ylim = c(.3, 4 + .3), 
     axes = F, xlab = NA, ylab = NA)       

abline(v = 0, lty = 2, col="grey") 

for (i in 1:length(dataret$dif)) {                                           
  points(dataret$dif[i], i, pch = 18, cex = 1.5)                             
  lines(c(dataret$low[i], dataret$up[i]), c(i, i))
  lines(c(dataret$low90[i], dataret$up90[i]),c(i,i), lwd=3)
  text(dataret$dif[i], i, dataret$V1[i], xpd = T, cex = 1, pos = 3)                    
}       

axis(side = 1)                                                                                         
mtext(side = 3, "(a) Support for Nuclearization", line = 0.9, cex=1.1)
box(col="grey")        


#Panel(b) Credibility of Extended Deterrence

plot(NULL,                             
     xlim = c(-80, 80),                     
     ylim = c(.3, 4 + .3), 
     axes = F, xlab = NA, ylab = NA)      

abline(v = 0, lty = 2, col="grey") 
                              
for (i in 1:length(datacet$dif)) {                                           
  points(datacet$dif[i], i, pch = 19, cex = 1.2)                             
  lines(c(datacet$low[i], datacet$up[i]), c(i, i))
  lines(c(datacet$low90[i], datacet$up90[i]),c(i,i), lwd=3)
  text(datacet$dif[i], i, datacet$V1[i], xpd = T, cex = 1, pos = 3) 
}


axis(side = 1)                                                                                        
mtext(side = 3, "(b) Credibility of Extended Deterrence ", line = 0.9, cex=1.1)

mtext(side = 1, "Changes in Support for Nuclearization \n(Threat-Declare)(%)", line = 3, at=-220, cex=1) 
mtext(side = 1, "Changes in Credibility of Extended Deterrence \n(Threat-Declare)(%)", line = 3, at=10, cex=1) 
box(col="grey")               


dev.off()






###################################################################################################
#### To Create Figure 3, calculate the difference in means and confidence intervals (Government Officials Survey)
####################################################################################################


##1) support for nuclearization 


mean(datae$support[datae$deploy==1])-mean(datae$support[datae$declaration==1]) ##-0.03
mean(datae$support[datae$threat==1])-mean(datae$support[datae$declaration==1]) ##0.19


##calculate confidence intervals 
set.seed(1000)
full.er <- bootTreat2(datae)

ed<-round(c(full.er$Deploy, quantile(full.er$boot[,1], c(0.025, 0.975))), digits=2)
ed#(-0.46, 0.40)
et<-round(c(full.er$Threat, quantile(full.er$boot[,2], c(0.025, 0.975))), digits=2)
et#(-0.26, 0.63)

ed90<-round(c(full.er$Deploy, quantile(full.er$boot[,1], c(0.05, 0.95))), digits=2)
ed90 #(-0.39, 0.33)
et90<-round(c(full.er$Threat, quantile(full.er$boot[,2], c(0.05, 0.95))), digits=2)
et90 # (-0.20, 0.57)





##2) credibility of extended deterrence 
mean(datae$credibility[datae$declaration==1]) #0.4545

mean(datae$credibility[datae$deploy==1])-mean(datae$credibility[datae$declaration==1]) ##-0.01
mean(datae$credibility[datae$threat==1])-mean(datae$credibility[datae$declaration==1]) ##-0.34


#calculate confidence intervals
set.seed(1000)
full.ec<-bootCredit2(datae)

ecd<-round(c(full.ec$Deploy, quantile(full.ec$boot[,1], c(0.025, 0.975))), digits=2)
ecd #(-0.46, 0.44)

ect<-round(c(full.ec$Threat, quantile(full.ec$boot[,2], c(0.025, 0.975))), digits=2)
ect #(-0.71, 0.04)

ecd90<-round(c(full.ec$Deploy, quantile(full.ec$boot[,1], c(0.05, 0.95))), digits=2)
ecd90 #(-0.39, 0.37)

ect90<-round(c(full.ec$Threat, quantile(full.ec$boot[,2], c(0.05, 0.95))), digits=2)
ect90 #(-0.64, -0.02)



##combine the results 




datae1<-as.data.frame(cbind(c("Threat", "Deploy")))
datae1$dif<-cbind(c(et[1], ed[1]))
datae1$low<-cbind(c(et[2], ed[2] ))
datae1$up<-cbind(c(et[3], ed[3] ))
datae1$low90<-cbind(c(et90[2],ed90[2]))
datae1$up90<-cbind(c(et90[3],ed90[3]))


datae1$dif<-datae1$dif*100
datae1$low<-datae1$low*100
datae1$up<-datae1$up*100
datae1$low90<-datae1$low90*100
datae1$up90<-datae1$up90*100


datae2<-as.data.frame(cbind(c("Threat", "Deploy")))
datae2$dif<-cbind(c(ect[1], ecd[1]))
datae2$low<-cbind(c(ect[2], ecd[2] ))
datae2$up<-cbind(c(ect[3], ecd[3] ))
datae2$low90<-cbind(c(ect90[2],ecd90[2]))
datae2$up90<-cbind(c(ect90[3],ecd90[3]))

datae2$dif<-datae2$dif*100
datae2$low<-datae2$low*100
datae2$up<-datae2$up*100
datae2$low90<-datae2$low90*100
datae2$up90<-datae2$up90*100





png("Figure3_FPA_elite.png", width=10, height=6,  units = 'in', res = 300)


m<-rbind(c(1,2))
layout(m)

## Panel (a) Support for Nuclearization 


plot(NULL,                            
     xlim = c(-80, 80),                    
     ylim = c(.3, 2 + .3), 
     axes = F, xlab = NA, ylab = NA)       

abline(v = 0, lty = 2, col="grey") 


for (i in 1:length(datae1$dif)) {                                          
  points(datae1$dif[i], i, pch = 18, cex = 1.5)                             
  lines(c(datae1$low[i], datae1$up[i]), c(i, i))
  lines(c(datae1$low90[i], datae1$up90[i]),c(i,i), lwd=3)
  text(datae1$dif[i], i, datae1$V1[i], xpd = T, cex = 1, pos = 3)                      
}       

axis(side = 1)                                                                                          # add bottom axis
mtext(side = 3, "(a) Support for Nuclearization", line = 0.9, cex=1.1)
box(col="grey")        

## Panel (b) Credibility of Extended Deterrence 

plot(NULL,                             
     xlim = c(-80, 80),                     
     ylim = c(.3, 2 + .3), 
     axes = F, xlab = NA, ylab = NA)       

abline(v = 0, lty = 2, col="grey") 

for (i in 1:length(datae2$dif)) {                                        
  points(datae2$dif[i], i, pch = 19, cex = 1.2)                              
  lines(c(datae2$low[i], datae2$up[i]), c(i, i))
  lines(c(datae2$low90[i], datae2$up90[i]),c(i,i), lwd=3)
  text(datae2$dif[i], i, datae2$V1[i], xpd = T, cex = 1, pos = 3) 
}


axis(side = 1)                                                                                        
mtext(side = 3, "(b) Credibility of Extended Deterrence ", line = 0.9, cex=1.1)

mtext(side = 1, "Changes in Support for Nuclearization \n(Treatment-Declare)(%)", line = 3, at=-220, cex=1) 
mtext(side = 1, "Changes in Credibility of Extended Deterrence \n(Treatment-Declare)(%)", line = 3, at=10, cex=1) 
box(col="grey")               

dev.off()






